def _periodic_bilinear_indices(
     grid: np.ndarray, values: np.ndarray, period: float
 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     n = grid.size
     if n < 2:
         raise ValueError("Interpolation grid must have at least 2 points")
 
     v = np.mod(values, period)
     i1 = np.searchsorted(grid, v, side="right")
     i0 = i1 - 1
     i0 = np.mod(i0, n)
     i1 = np.mod(i1, n)
 
     g0 = grid[i0]
     g1 = grid[i1]
 
     wrapped = i1 <= i0
     denom = np.where(wrapped, (g1 + period) - g0, g1 - g0)
+    numer = np.where(wrapped, np.mod(v - g0, period), v - g0)
     denom = np.where(denom == 0.0, 1.0, denom)
-    w = (v - g0) / denom
+    w = numer / denom
     w = np.clip(w, 0.0, 1.0)
     return i0, i1, w

 def load_coil_points(coil_file: str | Path) -> tuple[np.ndarray, np.ndarray]:
     """
     Load coil points from a SIMPLE/libneo-style coil file.
 
     Format:
         First line: integer N (number of points)
         Next N lines: x y z I
 
     Coordinates are expected to be in meters.
 
     Returns:
         points: (N, 3) array of xyz points [m]
         currents: (N,) array of currents
     """
     path = Path(coil_file)
     with path.open("r", encoding="utf-8") as f:
         header = f.readline().strip()
         if not header:
             raise ValueError(f"Empty coil file: {path}")
         try:
             n_points = int(header.split()[0])
-        except ValueError as exc:
-            raise ValueError(f"Invalid coil file header (expected N): {path}") from exc
+        except (ValueError, IndexError) as exc:
+            raise ValueError(f"Invalid coil file header (expected N points): {path}") from exc
 
         data = np.loadtxt(f, dtype=float, ndmin=2)
     if data.shape[0] != n_points or data.shape[1] < 4:
         raise ValueError(
             f"Coil file {path} expected {n_points} rows of x y z I, got {data.shape}"
         )
 
     points = data[:, 0:3].astype(float, copy=False)
     currents = data[:, 3].astype(float, copy=False)
     return points, currents

 def plot_heat_flux_with_coil_clearance(
     heat_map: WallHeatMap,
     *,
     coil_file: str | Path,
     chartmap_file: str | Path,
     min_clearance_m: float,
     ax=None,
     cmap: str = "hot",
     vmax: Optional[float] = None,
     show_colorbar: bool = True,
     ports: Optional[list[Port]] = None,
     forbidden_alpha: float = 0.35,
+    use_zero_current_separators: bool = True,
+    jump_factor: float = 25.0,
 ):
     """
 ... (clipped 10 lines)
     ax = plot_heat_flux_2d(
         heat_map,
         ax=ax,
         cmap=cmap,
         vmax=vmax,
         show_colorbar=show_colorbar,
         ports=ports,
     )
 
     wall = WallSurface.from_chartmap(chartmap_file)
-    coils = CoilSegments.from_file(coil_file)
+    coils = CoilSegments.from_file(
+        coil_file,
+        use_zero_current_separators=use_zero_current_separators,
+        jump_factor=jump_factor,
+    )
     clearance = clearance_map_on_heatmap_grid(
         wall,
         coils,
         heat_map.theta_centers,
         heat_map.zeta_centers,
     )
     forbidden = clearance < min_clearance_m
 
     forbidden_grid = forbidden.astype(float)
     ax.pcolormesh(
         np.degrees(heat_map.zeta_edges),
         np.degrees(heat_map.theta_edges),
-        forbidden_grid,
+        forbidden_grid.T,
         cmap="Greys",
         vmin=0.0,
         vmax=1.0,
         shading="flat",
         alpha=forbidden_alpha,
     )
     ax.set_title(
         f"Wall Heat Flux + Coil Clearance (d > {min_clearance_m:.3f} m)"
     )
     return ax

-for p, I in zip(points, currents, strict=True):
+for p, I in zip(points, currents):